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ABSTRACT 

Two  atmospheric  prediction  models,  based  upon  the  meteorological 
primitive  equations,  were  programmed  and  tested  on  the  CDC  6500  computer, 
A  frictionless ,  barotropic  model  was  integrated,  using  ten-minute  time  steps, 
at  500  MBS  for  periods  up  to  four  days .     Then,  an     n-level,  baroclinic  model 
for  an  inviscid,  adiabatic  atmosphere  was  integrated  for  brief  test  periods  at 
five  uniformly-spaced  zeta  (sigma)  surfaces  over  a  smooth  earth.    Both  models 
were  initialized  with  actual  data  fields  provided  by  FNWC  Monterey. 

Acceptability  was  based  on  measurements  of  energy  and  vorticity  para- 
meters, as  well  as  on  qualitative  assessments  of  output  fields.    The 
barotropic  500  MB  height  forecasts  were  found  to  be  comparable  to  the  FNWC 
(vorticity  model)  barotropic  forecasts.     Mean-square -vorticity  and  kinetic 
energy  were  suitably  conserved  for  integrations  up  to  four  days.    An  energy 
accumulation  in  the  smaller  range  of  scale  was  increasingly  evident  beyond 
two  days  into  the  forecast  period.     No  smoothing  was  performed  to  control 
this  accumulation,    A  limited  number  of  test  integrations  was  made  with 
the  five-level  model.    Computational  instabilities  were  observed  for 
forecasts  beyond  twelve  hours.    Interim  results  are  presented. 
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1.0  INTRODUCTION 

Disadvantages  arising  from  more  generalized  forms  of  the  vorticity, 
divergence,  and  thermodynamic  equations  have  overshadowed  the  advant- 
age of  relatively  large  integration  time  steps  in  filtered  prediction  models. 
At  the  same  time,  increased  knowledge  in  finite -difference  techniques, 
in  initialization  procedures,  and  in  boundary  condition  formulation,  have 
been  realized  in  a  decade  marked  by  unbelievable  increases  in  computer 
capacity  and  speeds.    These  trends  and  accomplishments  have  accelerated 
the  shift  to  the  so-called  "primitive  equations"  as  a  basis  for  atmospheric 
prediction  models,  both  operational  and  research. 

1.1  Finite  Differencing 

Whereas,  in  filtered  models,  an  integration  time  step  of  one  hour  is 
possible  without  risk  of  linear  computational  instability,  the  corresponding 
time  step  for  primitive -equation  integrations  is  approximately  ten  minutes. 
One  must  fulfill  the  von  Neumann  stability  criterion  in  order  to  avoid  this 
type  of  instability  in  linear  differential  equations.    A  second  type  of 
instability  may  arise  in  spite  of  careful  selection  of  time  step  and  grid 
mesh  length.    This  is  the  non-linear  computational  instability,  or 
"aliasing",  first  described  by  Phillips  (19  59).    Although  Phillips  treated 
non-linear  differential  equations,  Miyakoda  (1962)  showed  that  non-linear 
instability  can  also  occur  in  special  types  of  linear  equations. 

Aliasing  is  described  by  Lilly  (1965)  for  regions  in  which  some  arbi- 
trary continuous  function  is  being  represented  by  a  finite  number  of 
gridded  values.    It  is  shown  that  there  exists  a  maximum  and  minimum 
wavelength  which  can  be  resolved  by  any  grid  array,  and  that  waves 
outside  of  this  range  will  be  misrepresented  in  terms  of  the  permitted 
harmonics . 

13 


Phillips  controlled  this  phenomenon  by  smoothing.    It  may  be  reduced 
significantly,  however,  by  use  of  conservation  forms  of  the  difference 
analogs,  especially  those  which  preserve  the  average  wave  number.    The 
latter  prevent  the  cascade  of  energy  toward  smaller  ranges  of  scale.    The 
well-known  methods  of  Arakawa  (1966)  represent  those  which  conserve 
mean-square  vorticity,  total  kinetic  energy,  and  average  wave  number. 
Phase  errors  are  not  eliminated  by  conservation  forms  ,  but  their 
magnitudes  are  tolerable. 

Another  potential  difficulty  is  associated  with  the  phrasing  of  time 
derivatives,  which  have  taken  on  a  new  significance  with  the  advances 
made  in  space -differencing  methods.    The  work  of  Kurihara  (1965),  Lilly 
(1965),  and  Young  (1968),  is  illustrative  of  efforts  to  study,  and  possibly, 
to  reduce  the  "slow"  instability  associated  with  the  various  time -differenc- 
ing schemes.    The  "leapfrog"  scheme  of  Richtmyer  (1963)  has  been  used 
widely  in  spite  of  its  shortcomings.    These  may  be  listed  as  follows: 

a.  kinetic  energy  is  not  conserved, 

b.  during  extended  integrations,  non-linear  instability  may  develop 
in  the  absence  of  some  form  of  smoothing. 

c.  the  odd-  and  even-solution  sets  tend  to  proceed  independently  of 
each  other  because  of  the  alternating  sign  of  the  computational  mode  . 
(some  form  of  solution  averaging  may  be  helpful  to  eliminate  the  latter 
difficulty).    In  spite  of  these  shortcomings,  the  leapfrog  scheme  is  used 
because  of  the  implied  economies  and  its  general  accuracy. 

The  finite -differencing  scheme  ultimately  selected  for  a  given  model 
is  one  which  best  compromises  its  associated:  (1)  diffusive  character- 
istics; (2)  existence  and  properties  of  any  computational  mode;  (3) 
stability  properties;  and  (4),  the  computational  burden. 
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1.2      Boundary  Conditions 

The  boundary  conditions  for  any  model  must  be  tailored  both  to  the 
specific  equation  set  and  to  the  objective  of  the  integrations  (short  range 
forecast,  or  long-term  general  circulation  experiment).    The  primitive- 
equation  free -surface  barotropic  forecast  model  of  Shuman  and  Vanderman 
(1966)  is  indicative  of  the  effects  of  approximate  versus  correct  boundary 
conditions  on  the  forecast  fields.    Shuman  showed  that  the  RMS  velocity 
divergence  was  tolerably  conserved  for  periods  up  to  about  thirteen  days 
when  using  approximate  boundary  conditions,  but  that  the  model  "blew 
up"  rapidly  therafter.    By  using  the  correct  boundary  conditions  (which 
were  consistent  with  the  equations),  the  model  not  only  conserved  RMS 
velocity  divergence  from  the  outset,  but  also  prevented  further  deterioration 
of  the  forecast  which  used  the  approximate  conditions  for  the  first  thirteen 
days. 

The  important  point  is  that  incorrect  boundary  conditions  (as  well  as 
discretization  errors)  may  excite  gravity -inertial  waves  which,  in  the 
absence  of  some  dissipative  mechanism,  may  contaminate  the  forecast. 
See  Nitta  and  Hovermale  (1967)  for  an  interesting  discussion  of  these 
matters.    Matsuno  (1966)  reported  on  false  reflections  of  the  computation- 
al mode  from  an  outflow  boundary.    This  inward  propagation  of  errors  was 
especially  prominent  at  the  shorter  wavelengths.    The  rate  of  reflection 
decreased  with  increasing  wavelength  and  with  increasing  "order"  of 
continuity  at  the  lateral  boundaries.    Thus,  to  overcome  this  effect,  an 
increasing  number  of  interior  points  is  required  to  approximate  the 
boundary  value  of  an  advected  parameter  (if  a  suitable  assumption  or 
approximation  cannot  be  found  for  the  particular  problem). 
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1.3    Initialization  Procedures 

Initialization  procedures  are  probably  the  least  well-known  of  the 
requirements  for  model  development.    The  whole  problem  of  weather 
analysis  and  forecasting  is  concerned  with  the  mutual  adjustment  between 
the  mass  and  motion  fields  under  the  influence  of  rotation  and  gravity, 

One  very  interesting  result  arose  from  the  general  circulation  experi- 
ments of  Smagorinsky  (1965).    Patterns  of  precipitation  appeared  roughly 
the  same  regardless  of  the  initial  condition  of  vertical  motion,   provided 
that  other  initial  conditions  (such  as  temperature  and  the  non-divergent 
wind  component)  are  held  fixed.    Miyakoda  and  Moyer  (1968)  suggest 
that  for  periods  of  more  than  five  days  (for  the  integrations) ,  the  vertical 
motion  initialization  is  somewhat  arbitrary  because  of  the  deterministic 
nature  of  the  equations. 

Initialization  falls  into  two  classes:  (1)  extra pol at ive ,  and  (2) 
iterative  ,    In  either  case,  the  idea  is  to  reduce  the  noise  and  to  represent 
only  the  slowly  changing  meteorological  state,.    The  former  type  consists 
of  a  marching  integration  of  the  equations  to  arrive  at  this  "balanced" 
condition.    The  latter  type  stays  at  the  same  time  level,  but  iterates 
the  computations  until  the  equilibrium  state  is  reached,    Clearly,  the 
choice  depends  upon  the  model  purpose.     The  iterative  method  would 
have  an  advantage  for  users  requiring  accurate  short-term  forecasts  (of 
the  order  of  a  few  hours,  say),  whereas  the  extrapolative  method  would 
be  entirely  satisfactory  for  forecast  periods  beyond  the  time  required 
to  "settle  down"  (ten  hours  or  longer,  say). 

Thus  far,  primitive -equation  models  have  been  used  with  an  assort- 
ment of  initial  fields  -  some  real  and  some  analytical  constructions. 
Linearized  research  models  may  be  initialized  according  to  some  pre- 
selected analytical  distribution  with  a  known  solution 
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In  contrast,  Shuman  and  Hovermale  (1967)  chose  to  use  as  initial  data 
the  objective  analyses  of  geopotential  height  and  temperature  at  the  ten 
mandatory  pressure  levels  (from  1000  to  100  MBS).    Each  of  the  input 
analyses  was    expanded  from  the  JNWP  octagonal  grid  into  a  53  x  57 
rectangular  array  by  an  extra polative  scheme  which  consisted  of:  (1) 
finding  the  mean  value  of  the  analyzed  quantity  on  the  octagon  boundary, 
and  copying  this  value  into  the  non-overlapping  portion  (border  region) 
of  the  53  x  57  array;  and  (2)  "wedding"  the  two  regions  by  means  of  a 
"ring"  smoother,  while  holding  the  boundary  of  the  rectangular  array 
fixed.    By  this  process,  one  recognizes  that  the  equatorial  regions 
would  be  virtually  devoid  of  small  scale  components. 

The  removal,  or  reduction,  of  these  small-scale  features  is  clearly 
justified  for  several  reasons,  not  the  least  of  which  is  the  general 
inappropriate ness  of  the  so-called  balance  equation  in  the  tropics.    The 
difficulties,  moreover,  of  false  reflections  as  reported  by  Matsuno  (1966) 
would  be  minimized  in  such  "flat"  fields.    By  flattening  a  field  in  this 
way,  one  clearly  "increases  the  order  of  continuity"  at  the  expense  of 
accurate  specification  of  meteorological  detail  in  the  tropics  in  order 
to  preclude  possible  computational  difficulties. 

The  finite -difference  schemes,  boundary  conditions,  and  initializa- 
tion procedures  for  the  models  studied  in  this  paper  will  be  discussed 
in  detail  in  the  appropriate  sections  of  this  thesis. 

1.4    Coordinate  Systems 

Historically,  the  (x,y,p,t)  coordinate  system  was  shown  to  be 
superior  to  the  rectangular  system  in  many  respects,  but  it  suffered 
from  the  difficulties  arising  from  the  lower  boundary  condition, 
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Over  smooth  terrain,  one  generally  assumes  that  the  vertical  velocity 
vanishes  identically  at  the  lower  boundary;  whereas,  over  rough  terrain, 
one  assumes  that 

u>±  \V.  V-rr 

(or  some  equivalent  approximation)  where  IT  is  the  terrain  pressure - 
In  Phillip's  (x,y,<T,t)  sigma  system,  the  particular  problem  just 
mentioned  is  overcome,  but  the  sigma  system  is  not  entirely  free  of 
practical  difficulties.     Kurihara  (1968)  reported,  for  example,  on  problems 
arising  from  (large)  terrain-pressure  gradients . 

In  this  research  project,  the  barotropic  model  was  written  in 
pressure  coordinates,  but  the  baroclinic  model  was  based  upon  a 
generalized  "zeta"  system  which  was  proposed  by  R.T.  Williams 
(personal  communication).    In  this  modified  sigma  system,  the  vertical 
coordinate  zeta  is  related  to  Phillip's  sigma  by  an  arbitrary  functional 
relationship 


«"=■  +  (S") 


where,  it  is  clear,  one  may  achieve  the  desired  "linearity"  in  the 
vertical  parameter  while  allowing  for  arbitrary  placement  of  the 
computational  surfaces. 
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2  .  0      THE  BAROTROPIC  EXPERIMENT 
2.1     Introduction 

The  matter  of  testing  a  simple  barotropic  model  prior  to  engaging  in 
more  complex  simulations  is  a  desirable  approach.    The  work  of  Houghton 
et  al  (1966)  and  Shuman  and  Vanderman  (1966)  is  cited  in  this  regard. 

Houghton  conducted  experiments  with  various  finite -difference 
formulations  for  the  Coriolis  term  and  boundary  conditions  in  order  to 
observe  sensitivities  in  long-term  solutions.    Solution  stability  was 
achieved  without  a  large  amount  of  smoothing,  but  some  was  necessary 
in  order  to  control  truncation  errors.    Charney  (1965)  had  suggested  an 
upper  limit  (several  weeks)  for  long-term  integrations  because  of 
inaccuracies  in  initial  conditions.    Houghton  allowed  that  stability  and 
truncation  effects  may  be  of  equal  importance  to  the  validity  of  long- 
term  solutions  as  to  the  growth  of  initial  error. 

In  this  project,  the  author  has  speculated  that  the  small  errors,  or 
inconsistencies,  in  the  initial  conditions,  as  well  as  the  effects  of 
truncation  and  discretication,  will  not  be  of  major  concern  for  short- 
term  integrations.    There  does  appear  to  be  some  basis  for  more  concern 
about  the  former  in  one-  and  two-day  forecasts. 

Shuman  and  Vanderman  (1966)  also  tested  a  primitive -equation  free- 
surface  barotropic  model  in  a  tropical  channel  experiment.    They 
concluded  that  approximate  boundary  conditions  (which  were  inconsistent 
with  the  difference  analogs)  were  entirely  satisfactory  for  integrations 
up  to  about  twelve  days.    The  philosophy  behind  Shuman's  choice  was 
that  the  conditions  provided  for  reflection  of  pure  gravity  waves  at  the 
walls.    The  inertial    effects  (due  to  the  presence  of  "f"),  it  was  reasoned, 
could  be  ignored  for  a  limited  time. 
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The  finite-difference  analog  referred  to  by  Shuman  as  the  "semi- 
momentum"  form  apparently  contains  a  large  amount  of  smoothing 
because  of  the  generous  use  of  the  superbars.    This  may  account  not 
only  for  the  success  of  the  long-term  integrations  without  explicit 
smoothing,  but  also  for  the  degree  of  concern  about  "unsmoothing" 
before  use  of  the  balance  equation. 

2 . 2    Purpose  of  the  Barotropic  Experiment 

The  first  phase  of  thesis  research  consisted  of  the  derivation, 
programming  (CDC  6500),  and  evaluation  of  a  barotropic  free-surface 
model.    This  was  test  run  at  the  500  MB  level  using,  as  initial  data, 
the  actual  objectively-analyzed  hemispheric  data  fields  provided  by 

FNWC  Monterey. 

The  author,  after  noting  the.  many  achievements  of  investigators 
cited  earlier,  felt  that  the  obstacles  remaining  for  reasonably-accurate 
short-term  integrations  of  real  data  fields  were  more  practical  than 
theoretical.    If  would  instructive,  that  is,  to  see  if  approximations 
and  assumptions  found  suitable  in  analytical  test  integrations  were 
equally  suitable  for  integrations  with  real  data  fields, 

The  purpose  of  the  barotropic  experiment  was  to  note  the  suitability 
of  certain  boundary  conditions,  to  determine  appropriate  initialization 
procedures,  and  to  test  the  modified-Arakawa  finite -difference  analogs 
(as  proposed  by  R.T.  Williams)  for  a  model  which  might  be  expected  to 
reasonably  simulate  atmospheric  flow  patterns  at  500  MBS  for  periods 
up  to  four  days. 
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2.3     The  System  of  Prognostic  Equations 

The  two-dimensional  differential  equations  of  momentum  and 
continuity  for  a  free -surface  barotropic,  inviscid  atmosphere  may  be 
approximated  on  a  polar  stereographic  projection  as  follows: 

a.  momentum  equation  in  east-west  direction: 

a-b       L         3  ax       v    ax        au  ;i 

-  rm       r-     —  \    4-  3-  (  viii—  \  L  2.3.1 

L<)x^   rm  )        a^\    ^  )J 

b.  momentum  equation  in  north-south  direction: 

-W$  |-(ilid£\  x-iL/  Uir  xt  (2.3.2) 


c.    continuity  equation: 


(2.3.3) 


The  reader  is  referred  to  Appendix  A  for  the  development  of  this  set 
of  equations  in  flux  form. 

The  computational  boundary  was  taken  to  be  mid-way  between  the 
outer  ring  of  grid  points  and  the  first  interior  ring  of  points .    This 
permitted  the  practical  application  of  boundary  conditions  required  for 
conserving  the  square  of  any  parameter  being  advected,  while  at  the 
same  time,  permitting  non-zero  normal  wind  components  at  the  first 
interior  ring. 
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Geostrophic  flow  was  assumed  for  the  normal  wind  component  along 
the  first  interior  ring  of  points  because  of  the  computational  convenience 
By  so  doing,  the  problem  of  computing  the  correct  values  of  "h"  along 
the  outermost  ring  is  avoided. 

The  advective  terms  led  to  the  following  boundary  conditions .    At 
the  computational  boundaries  parallel  to  the  y-axis, 

— n* 


o 


(2.3.4) 


while,  at  the  boundaries  parallel  to  the  x-axis , 

— r Y 


(2.3.5) 


For  convenience,  the  tangential  wind  component  at  the  first  interior 
ring  was  copied  into  the  adjacent  grid  points  on  the  outermost  ring.    In 
a  similar  manner,  the  values  of  "h"  were  also  shifted  "out"  for 
convenience. 

It  will  be  seen  that  since  the  minimum  value  of  the  Coriolis 
parameter  permitted  was  that  value  corresponding  to  "f"  at  15     latitude, 
the  map  factor  advection  terms  (which  are  small  anyway)  vanish  south 
of  that  latitude.    The  local  change  of  "utv*  and  "vh"  depended,  therefore, 
on  the  horizontal  flux  divergence  terms  at  the  first  interior  ring  (since 
geostrophic  flow  was  assumed).    Single  (space)  differences  were  used 
as  appropriate  in  the  continuity  equation  along  the  first  interior  ring. 
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Initialization  of  wind  components  consisted  of  solving  the  so-called 
linear  balance  equation  of  the  form 


V  <§  -    7  •  (+^7+) 


(2.3.6) 


using  an  FNWC  500  MB  height  field.    The  "first  guess"  streamfunction 
values  were  taken  to  be 

T  -T  (2.3.7) 

T 

o 
where    h    was  the  initial  height  value  at  each  point. 

The  integration  time  step  was  ten  minutes.    Subsequent  to  an 

initial  forward  (single)  step,  a  conventional  centered  difference  was 

used. 

2.4     Finite  Difference  Forms 

A  modified-Arakawa  finite -difference  scheme  was  employed.    In 
this  scheme,  the  momentum  and  continuity  equation  analogs  were 
expressed  in  the  forms  shown  below: 

a.    momentum  equation  in  east-west  direction: 

M-M  M-l  (        r  r         \ 


*S,I(VS^^ 


^^K.^r^^.^ji  <2-^ 
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b.     momentum  equation  in  north-south  direction: 


M 


^jKj+.-W..,Nni-  (2-4.2) 


J 

where,  for  any  advected  parameter,  [b  ,  we  have  defined  the  operator 


(%:ti)M,J% 


(2.4.3) 


c.    continuity  equation 


where  "n"  denotes  the  time  level. 
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In  an  attempt  to  reduce  some  of  the  small  scale  irregularities  of  the 
initial  height  analyses  near  the  computational  boundary,  a  five -point 
smoother  routine  was  written  to  produce  the  following  effect.    First, a 
mean  height  value  (say)  was  calculated  for  the  outermost  grid  ring. 
This  value  was  set  at  all  of  the  outermost  ring  points.    The  smoothing 
coefficient  was  latitude  dependent  in  the  sense  that  the  smoothing 
increased  from  zero  at  20    latitude  to  a  degree  sufficient  to  remove 
all  perturbations  at  the  outermost  ring.    This  procedure  was  deemed  to 
be  consistent  with,  and  justified  by,  the  linear  balance  assumption 
in  the  tropics. 

2.5     Research  Procedures  and  Results 

The  barotropic  model  was  tested  on  actual  FNWC  500  MB  height 
analyses  for  the  Northern  Hemisphere,  rather  than  on  some  idealized 
distribution  for  which  some  analytical  solution  was  known.    This  less 
rigorous  test  was  consistent  with  the  purpose  of  the  study  -  to  see  if 
some  of  the  approximations  found  suitable  for  analytical  tests  were 
equally  suitable  for  use  with  actual  data  fields  (and  their  much  larger 
frequency  ranges). 

The  judgement  as  to  solution  stability  or  acceptability  was  based 
on  certain  diagnostic  parameters,  such  as  the  square -vorticity 
parameter  (SVP)  or  mean  kinetic  energy  (KE) ,  as  well  as  on  the  meteor- 
ological appearance  of  the  output  fields . 

The  first  test  consisted  of  computing  all  parameters  and  terms  as 
they  appeared  in  the  equations  (without  space  averaging).    Then,  the 
model  was  rerun  using  space  averaging  in  a  manner  similar  to  that 
used  by  Shuman  (1966).    It  was  determined  that  the  SVP  grew  at  an 
intolerably  large  rate. 
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When  no  space  averaging  was  performed,  the  SVP  growth  rate  was  an 
order  of  magnitude  smaller,,    In  spite  of  the  lesser  rate,  however,  energy 
did  tend  to  accumulate  in  the  smaller  range  of  scale,  particularly  in 
tropical  regions . 

A  second  test  consisted  of  measuring  the  SVP  for  a  run  in  which  the 
integration  cycle  was  restarted  every  72  time  steps,  but  with  no  space 
averaging.  This  was  an  attempt  to  measure  the  effects  of  non-linear 
interactions  between  the  odd-  and  even-solution  sets.  The  result  was 
that  the  SVP  grew  at  approximately  the  same  rate  for  both  runs,  but  the 
SVP  took  a  small  step  upward  in  value  immediately  following  the 
restart. 

Figure  1  depicts  the  aforementioned  SVP  growth  rates.    The  upper 
curve  (A)  shows  the  rate  for  the  case  when  space  averaging  was  used, 
but  no  restarting  was  performed.    The  lower  curve  (B)  shows  the  SVP 
rate  for  both  cases  in  which  no  space  averaging  was  performed  (no 
attempt  was  made  to  indicate  the  small  upward  jogs  after  each  restart, 
since  they  were  negligible  on  this  scale).    The  larger  rate  (A)  was 
judged  to  be  the  result  of  some  rather  serious  imbalances  caused  by 
not  "unsmoothing"  initial  fields  in  the  manner  reported  by  Shuman  (1966, 
1967).     The  much  smaller  rate  (B)  was  most  likely  caused  by  small 
imbalances  arising  from  the  use  of  the  linear  balance  equation,  and  by 
the  gradual  accumulation  (which  could  be  either  a  real  or  computationally- 
induced  cascade)  of  energy    in  the  small  scale  features.    No  smoothing 
was  performed  in  case  B, 
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Figure  2.    Difference  Field  (in  meters)  between  the  odd-  and  even- 
s  Nation  sets  after  144  ten-minute  time  steps  in  barotropic  experi- 
ment.   This  area  shown  represents  a  lOd  x  lOd  subset  of  the  total 
:.?mispheric  field.     It  is  idealized  to  portray  the  typical  scale 
-iid  magnitude  of  such  a  difference  field.    In  general,  the  range  of 
;c~>\?  of  features  is  from  2d  tc<  4d,  with  magnitudes  generally  less 
r.r  an  five  meters  . 
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A  third  test  was  made  to  measure  the  "departure"  of  the  two 
solution  sets.    Figure  2  illustrates  the  (idealized)  scale  and  magnitudes 
of  the  difference  field  after  144  time  steps.    Each  side  of  the  figure  is 
ten  mesh  lengths  long  (about  3810  kilometers),  and  is  representative  of 
the  entire  hemispheric  difference  pattern.    In  general,  the  magnitudes 
of  the  centers  are  five  meters  or  less  after  twenty -four  hours  of  forecast 
period.    This  may  be  compared  to  the  magnitude  of  the  500  MB  diurnal 
height  oscillation,  which  is  (roughly)  an  order  of  magnitude  larger  than 
this  difference  field. 

The  mean  kinetic  energy  was  computed  for  each  time  step.    No 
attempt  was  made  to  plot  these  distributions  because  the  variation  was 
negligible  in  the  cases  in  which  no  space  averaging  was  performed. 

One  test  was  attempted  in  which  the  latest  odd-  and  even-step 
wind  component  and  height  fields  were  averaged  (as  appropriate)  after 
72  time  steps,  and  then  restarted.    The  calculations,  however,  blew 
up  immediately  after  restarting  the  cycle.    Thus,  it  would  appear  as 
though  the  "average"  wind  components  and  "average"  heights 
introduce  a  considerable  trauma  in  the  system,  in  general.    It  should 
be  noted  that,  at  all  other  times,  the  restarting  procedure  involved 
using  the  latest  fields  rather  than  averaged  fields. 

Charts  A  through  E  depict  the  initial  500  MB  analysis  and  daily 
forecast  fields  out  to  four  days,  respectively.  The  model  tends  to 
overdevelop  both  highs  and  lows ,  but  it  is  felt  that  the  phenomenon 
could  be  easily  controlled.  Note  also  that  the  amplitudes  of  small 
scale  features  increase  perceptibly  (with  time),  particularly  in  low 
latitudes.    No  attempt  was  made  to  control  this  accumulation  of  energy. 
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The  barotropic  primitive -equation  24-hour  500  MB  forecast  fields 
were  compared  with  the  barotropic  filtered  model  output  fields  provided 
by  FNWC.    The  height  error  fields  revealed  that  the  PE  model  output 
compared  favorably  with  the  operational  FNWC  model  in  the  case 
tested.    Charts  F  and  G  illustrate  the  height  error  patterns  for  the 
PE  and  filtered  models,  respectively,  for  the  24-hour  period  ending 
0000Z,   25  April  1968.    It  was  beyond  the  intent  of  this  research  to  make 
more  than  a  perfunctory  comparison  of  these  two  models. 

In  summary,  it  can  be  said  that  the  initialization  through  use  of  the 
linear  balance  equation  was  satisfactory  when  small-scale  features 
were  removed  from  the  boundary  regions  by  a  latitudinally-varying 
filter.    The  boundary  conditions  and  finite -difference  forms  used  in 
the  experiment  were  both  economical  and  stable  for  integrations  up  to 
four  days.    Although  the  smaller-scale  features  tended  to  increase  in 
amplitude,  it  was  felt  that  the  phenomenon  could  be  easily  controlled. 
The  results  would  justify  further  research  on  a  barotropic  model 
initialized  with  real  hemispheric  data  fields. 
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CHART  A 


Chart  A.      Initial  500  MB  Height  Analysis  for  0000Z,   24  April  1968. 
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CHART  Bl 


Chart  B.      P.E.  Barotropic  500  MB  24-Hour  Height  Forecast  from 
0000Z,    24  April  1968. 
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Chart  C.      P.E.  Barotropic  500  MB  48-Hour  Height  Forecast  from 
0000Z,   24  April  1968. 
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CHART  D 


Chart  D.      P.E.  Barotropic  500  MB  72-Hour  Height  Forecast  from 
0000Z,   24  April  1968. 
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Chart  E.      P.E.  Barotropic  500  MB  96-Hour  Height  Forecast  from 
0000Z,   24  April  1968. 
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CHART  F 


Chart  F.      P.E.  Barotropic  24-Hour  500  MB  Height  Error  Field 
(in  meters)  for  period  ending  0000Z,   25  April  1968. 
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CHART  G 


Chart  G.      FNWC  (Vorticity  Model)  Barotropic  24-Hour  500  MB 
Error  Field  (in  meters)  for  period  ending  0000Z, 
25  April  1968. 
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3.0  THE  BAROCLINIC  EXPERIMENT 

3.1  Introduction 

A  generalized  Fortran  program  was  written  for  integrating  an  n-level 
baroclinic  primitive -equation  set  for  a  frictionless ,  adiabatic  atmos- 
phere.   Because  of  its  generality,  much  of  the  possible  computational 
economy  was  sacrificed.    As  a  result,  the  frequent  use  of  auxiliary 
memory  (disc  file)  caused  the  expenditure  of  about  twice  as  much 
peripheral  processor  time  as  central  processor  time  per  time  step . 

To  check  out  and  test  run  the  model,  the  input-output  section 
was  written  for  five  uniformly-spaced  zeta  (sigma)  surfaces  only,    The 
vertical  coordinate  (zeta)  may  be  given  by  the  transformation 

<T=     -f  (S)  (3.1,1) 

where    (T     is  defined  as  the  ratio  of  pressure  to  that  at.  the  lower 
boundary  (1T);  that  is, 

T  (3.1,2) 

Differentiation  with  respect  to  sigma  leads  to 

*"   \»)  (3-K3) 

which,  in  turn,   suggested  the  definition  of  the  so-called  "vertical 
map  factor" 


Sm-rt£T' 


V<>5"/  »•!.« 
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This  leads  to  the  derivative  transformation 

5?"    sa?  (3-U5) 

where  A  is  an  arbitrary  dependent  variable.    In  the  special  case  used 
for  evaluating  the  model,  the  transformation  was  linear  and  the  vertical 
map  factor,  s,  was  equal  to  unity.    Moreover,  since  zeta  increased 
upward,  the  vertical  velocity  was  defined  positive  upward  according 
to  the  relation 

ur  =   -  <r 

(3.1.6) 

It  is  well-known  that  Phillip's  sigma  system  (and  therefore,  the 
zeta  system)  has  both  advantages  and  disadvantages.    Its  main 
advantages  are: 

a.  the  absence  of  one  of  the  chief  difficulties  in  pressure 
coordinates  -  that  of  needing  to  define  the  horizontal  pressure  gradient 
at  sea  level.      For  models  that  include  terrain,  the  generation  of 
ficticious  pressure  gradients  is  avoided  because  the  terrain  pressure 

is  observed  rather  than  calculated  through  some  reduction  scheme. 

b.  the  function  representing  vertical  motion  vanishes  identically 
at  all  material  surfaces. 

Two  of  the  disadvantages  of  the  sigma  system  are: 

a.  that  observational  data  are  generally  taken  at  pressure 
surfaces,    Thus,  sigma  surface  parameters  must  be  generated  by  some 
interpolation/extrapolation  scheme . 

b.  that  computational  difficulties  may  arise  because  of  large 
space  variations  of  terrain  pressure  in  certain  regions  (see  Kurihara, 
1968). 
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3.2  Purpose  of  the  Baroclinic  Experiment 

Realizing  that  time  would  not  permit  the  inclusion  of  diabatic  heating 
and  diffusion  terms  in  the  equations,    the  goal  was  made  to  program 
an  n-level  model  for  an  adiabatic,  inviscid  atmosphere  on  a  smooth 
earth;  and  to  test  the  model  at  five  uniformly-spaced  zeta  surfaces  (in 
the  time  remaining) .    Provisions  were  made  in  the  Fortran  program  to 
incorporate  the  remaining  terms  later  without  need  for  any  major 
re  programming  effort,, 

As  in  the  barotropic  experiment,  it  remained  to  continue  the 
evaluation  of  boundary  conditions  and  initialization  procedures,    The 
modified-Arakawa  conservation-type  difference  scheme  had  to  be  tested 
in  the  multi-level  version.    In  addition,  the  author  was  aware  of  the 
buildup  of  energy  in  small-scale  features  in  the  barotropic  model,  and 
desired  to  observe  the    manifestations  of  frictionless  flow  in  the  more 
complex  model. 

3.3  The  System  of  Prognostic  Equations 

The  equation  set  for  an  inviscid,  adiabatic  atmosphere  may  be 
approximated  on  a  polar  stereographic  map  projection  as  follows: 
a,    east-west  momentum  equation: 
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b.    north-south  momentum  equation: 


ainr=  -  ^  \  #(™iT\ + 2.  ( ££  V 


_^S  3(uJur) 


c.    continuity  equation: 


d.    hydrostatic  equation: 


5?  ~    s<r 


e.    thermodynamic  equation: 


[axV  mijTajVvn' 


** 


Tps  a(ur  i ) 

as" 


(3.3.2) 


(3.3.3) 


(3.3.4) 


(3.3.5) 


where  all  symbols  have  been  defined  previously. 

Figure  3  contains  a  diagram  of  computational  levels  with  notations 
as  to  which  variables  are  predicted  or  computed  at  those  levels. 
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Level 

Computed  Variables        (Notational 
Zeta  (at  each  level)  Subscript)  Sigma 

1.0    _ _ ___^_^___^___    0.0 

0.9 J±LV±L>* 5 0.1 


0.8 - 2 0.2 


w  3 

0.6 ■    0.4 


0.5 1J-' 0.5 

w                                           2 
0.4 0.6 

0.3 'dlV^'i- 2 0.7 

w                                         1 
0.2 0.8 

u, v,T, z                                    1 
0.1 0.9 

_    n  Boundary  Pressure  (pi) 


Figure  3.     Diagram  it  Zeta  (Sigma)  Levels.     Notations  are  made  to 
indicate  the  subscripts  used  for  the  level,  as  well  as  to  show  the 
variables  predicted  or  computed  at  each  level.     Note  that  "z"  has 
been  used  interchanaeably  with  phi  (geopotential) ,  and  that  "pi" 
ha-;  been  used  for  the  1  iwer  boundary  pressure. 


42 


3.4     Finite-Difference  Forms 

The  difference  equations  corresponding  to  (3„3)  are  based  on  the 
Arakawa  technique.  It  will  be  noted  that  the  set  (3.3)  are  similar  in 
form  to  those  used  earlier  by  Smagorinsky  (1965). 

The  difference  scheme  avoids  non-linear  instability  by  requiring 
that  the  advective  terms  conserve  the  square  of  the  advected  property 
(when  summed  over  the  entire  grid),  assuming  continuous  time 
derivatives.    The  total  energy  is  conserved  because  of  requirements 
placed  upon  the  vertical  differencing;  specifically,  upon  the  special 
form  of  the  hydrostatic  equation.    Assuming  that  the  zeta  surfaces  are 
uniformly-spaced,  the  lowest  phi  distribution  is  given  by 

RA5  /    T 


(3. 4. la) 


and  for  subsequent  levels,  one  uses 


(3. 4. lb) 


It  should  be  noted  that  a  more  complicated  form  of  the  hydrostatic 
equation  is  required  if  the  surfaces  are  not  uniformly-spaced.    The 
reader  is  referred  to  Haltiner  (3)  for  a  complete  discussion  of  the 
conservation  properties  of  this  difference  scheme. 
The  east-west  momentum  equation  is  given  by 

(A 

where  (*)  denotes  the  point  (i,j,k), 
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(3.4.2) 


Similarly,  the  north-south  momentum  equation  is  written 


-<' 


The  lower  boundary  pressure  is  computed  from 

M-l 


/v\ 


(3.4.3) 


-tr 


(3.4.4) 


J 


J 


J 


where,  for  the  entire  column,  one  evaluates  the  local  change  as 


where,  for  notational  convenience,  one  defines  the  quantities 


(3.4.5) 


ot  = 


/Yv\ 


and 


/3»  S 

l-*  AAA 


The  thermodynamic  equation  takes  on  the  form 


wt;  =  (1ff)  +  a  tt  - <cw  ^*  [-ii  C^+Wfc  )i ; 


^%*\^%r^Y^r*A}} 


/A 


(3.4,6) 
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Note  that  for  any  advected  parameter  an  operator  has  been  defined 

3. 


♦flfc 


^W-^A+t.) 


ic-r  fc       k-i 


(3.4.7) 


li 


Finally,  vertical  motions  for  non-material  surfaces  are  obtained 
from  a  form  of  the  continuity  equation 

3.5     The  Computational  Sequence 

For  the  initialization  portion  of  the  program,  the  sequence  of  steps 
is  as  follows: 

a.  compute  map  factor  and  Coriolis  parameter. 

b.  read  input  fields  from  magnetic  tape: 

(1)  sea  level  pressure. 

(2)  temperatures  at  850,  700,   500,  400,  300,  200,  and 
100  MBS. 

c.  convert  temperatures. 

d.  interpolate  p-surface  temperatures  to  build  zeta-surface 
temperatures,  and  smooth  fields  in  tropics. 
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e .  form  T-pi  fields  , 

f.  integrate  hydrostatically  to  obtain  phi  fields. 

g.  solve  linear  balance  equation  at  each  level. 
h,    compute  u-  and  v-components  for  each  level, 
i.    form  u-pi  and  v-pi  fields. 

For  the  steps  involved  in  each  time  step  after  initialization  is 
completed,  the  following  are  listed: 

a.  compute  local  change  in  lower  boundary  pressure. 

b.  compute  vertical  motion  fields. 

c.  compute  next  lower  boundary  pressure, 

d.  compute  forcing  functions  for  momentum  and  thermodynamic 
equations  for  this  level. 

e.  integrate  to  obtain  next  u-pi,  v-pi,  and  T-pi  fields  for 
this  level. 

f.  form  next  u,  v,  and  T  fields  for  this  level, 

g.  repeat  steps  d.  through  f.  for  remaining  levels, 
h.    compute  next  phi  fields. 

i,     (optional)  compute  diagnostic  parameters » 

It  should  be  noted  that  the  integration  cycle  was  restarted  with  a 
single  forward  time  step  every  six  hours  in  the  baroclinic  experiment, 

3.6    Boundary  Conditions 

The  boundary  placement  is  identical  to  that  used  in  the  barotropic 
experiment.  As  in  the  barotropic  experiment,  the  advective  terms  led 
to  the  following  boundary  conditions . 
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At  the  computational  boundaries  parallel  to  the  y-axis, 

-X 

UTT  \     _ 

o 


(3.6.1a) 


and 


^r-i    =      TT  ~&  (3.6.1b) 


while,  at  the  boundaries  parallel  to  the  x-axis, 

Y 


l>£TL\_-a  (3.6.2a) 

VYV1 ) 


and 

5$  ^>TT 


-   o 


^     ^ 


(3.6.2b) 


The  practical  imposition  of  boundary  conditions  was  varied  slightly 
from  one  test  to  another.    For  the  most  part,  however,  geostrophic 
flow  was  assumed  at  the  first  interior  ring  in  order  to  avoid  the 
necessity  of  estimating  the  values  of  "phi"  and  "pi"  along  the  outer- 
most ring.    One-sided  space  derivitives  were  used  in  the  continuity 
and  thermodynamic  equations  (as  appropriate)  for  computations  along 
the  first  interior  ring. 

At  the  earth's  surface  and  at  the  top  of  the  atmosphere,  the  vertical 
velocity,  w,  vanished  identically. 
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3.7     Initialization  Procedures 

The  input  fields  for  the  baroclinic  experiment  were  actual  data 
fields  as  analyzed  on  the  FNWC  63  x  63  hemispheric  grid.    To  remove 
most  of  the  disturbance  components  from  the  boundary  regions,  a 
latitudinally-de pendent  smoother  of  the  form 

was  employed  to  compute  the  coefficient.    The  coefficient  vanished 
north  of  30     latitude, 

Initialization  of  the  geopotential  and  wind  fields  was  accomplished 
by  use  of  the  linear  balance  equation  for  a  sigma(zeta)  surface,     (It 
should  be  noted  that  the  terms  arising  from  the  presence  of  the 
sphericity  terms  in  the  momentum  equation  were  not  retained.)    The 
balance  equation  took  on  the  form 

V>  =   H^  +  RT  j  VW  -R^W^f ) 


4 


"T  J  (3.7,2) 


as  a  result  of  the  divergence  of  the  two-termed  pressure  gradient  force 
in  sigma  coordinates. 

Appendix  C  contains  a  derivation  of  the  divergence  equation  in 
sigma  coordinates. 
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3.8      Research  Procedures  and  Results 

The  generality  of  the  Fortran  program  led  to  a  very  inefficient  use 
of  the  CDC  6500  computer  at  the  Fleet  Numerical  Weather  Central.    For 
the  five -level  version  tested,  approximately  180  disc  accesses  per  time 
step  were  required  (for  storage  of  the  history  variables,  and  temporary 
storage  of  intermediate  results).    Because  of  the  extremely  heavy 
burden  of  high-priority  operational  programs  being  run  on  the  FNWC 
computer  in  the  period,  October  1968  to  December  1968,  this  phase 
of  research  was  limited.    In  addition  to  "debugging"  the  2000-instruction 
program,  a  few  integrations  were  run  out  to  twelve  hours  in  order  to 
make  a  few  assessments  of  the  output  fields  and  measurements  of  the 
conservation  parameters . 

Interim  results  of  the  baroclinic  experiment  indicate  the  generation 
of  static  instabilities  in  the  uppermost  layers,  and  a  subsequent  rapid 
deterioration  in  the  fields .     No  attempt  was  made  to  constrain  the 
lapse  rates  of  predicted  temperatures. 

The  output  fields  at  the  sixth  hour  into  the  forecast  did  not  show 
any  indication  of  computational  deterioration. 

In  almost  every  instance,  the  computational  blow-up  was 
accompanied  by  a  catastrophic  potential-to-kinetic  energy  conversion 
in  the  uppermost  layers . 

A  continuing  series  of  test  runs  are  being  made  by  the  author  in 
an  effort  to  ascertain  the  source  of  the  computational  difficulties. 
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4.0  CONCLUSIONS  AND  RECOMMENDATIONS 

4.1  The  Barotropic  Experiment 

The  barotropic  experiment  was  intended  to  test  the  suitability  of 
boundary  conditions,  initialization  procedures,  and  the    phrasing  of 
the  finite -difference  analogs.    In  addition,   some  information  as  to 
the  degree  of  smoothing  of  the  initial  data  fields  in  the  tropical  regions 
was  sought. 

Initialization  through  use  of  the  linear  balance  equation  was 
satisfactory  when  most  of  the  disturbance  component  was  removed 
at  the  boundary  region  through  use  of  a  latitudinally-varying  smoother 
routine . 

The  boundary  conditions  and  finite -difference  forms  were  found 
to  be  both  economical  and  stable  for  integrations  up  to  four  days.   (No 
runs  were  made  beyond  four  days.) 

Although  energy  tended  to  accumulate   in  the  smaller  ranges  of 
scale,  it  was  felt  that  this  phenomenon  could  be  controlled. 

The  results  of  this  investigation  would  justify  further  research  on 
a  barotropic  primitive -equation  model  that  was  initialized  with  real 
hemispheric  data  fields.    The  conservation  properties  of  the  difference 
scheme  place  the  model  at  an  advantage  over  the  operational  FNWC 
model  for  intermediate  length  integrations  (three  to  seven  days). 

4.2  The  Baroclinic  Experiment 

Since  most  of  the  test  runs  made  indicated  that  instabilities 
developed  in  the  first  eight  to  thirteen  hours  of  the  forecast  period, 
more  work  is  required  to  ascertain  the  nature  and  cause (s)  of  the 
computational  instabilities. 
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Preliminary  results  would  seem  to  indicate  that  the  vertical  velocity 
fields  were  not  satisfactory,  generally.    This  was  felt  to  be  particularly 
true  in  the  upper  atmosphere,  where  sudden  transformations  of  potential 
to  kinetic  energy  were  linked  to  unrealistic  vertical  transports,  and  the 
ensuing  static  instabilities. 

The  generalized  Fortran  program  should  be  useful  as  a  vehicle  for 
future  graduate  research  work  at  the  Naval  Postgraduate  School. 
Provisions  have  been  made  for  arbitrarily  varying  the  number  and 
location  of  zeta  levels.      In  addition,  heating  and  friction  terms  may 
be  added  without  need  for  extensive  re  programming. 

An  extended  core  system  (ECS)  was  added  to  the  CDC  6500  computer 
in  late  December  1968,  and  will  provided  an  additional  500,000  words 
of  fast-access  auxiliary  memory.    With  little  additional  programming, 
this  should  greatly  reduce  the  "field  handling"  time,  and  permit  longer 
test  runs . 
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APPENDIX  A 
Flux  Forms  of  Barotropic  Equations  in  Pressure  Coordinates 

The  flux  forms  of  the  barotropic  primitive -equation  set  in  pressure 
coordinates  may  be  "transformed"  from  the  standard  forms  for  a  polar 
stereographic  projection.    These  forms  to  be  transformed  are  given  by: 


(A.  la) 


(A.  lb) 


(A.  2) 


The  procedure  consists  of  showing  that  local  change  terms  may  be 
written  as: 


(A.  3a) 


(A .  3b) 


Similarly,  the  advective  terms  lead  to  the  forms, 


h.  =     C21  f  \  3  /  box  +  S>j  t,Mirx| 


m^  I  ^4-  **       ^  'J 


/ 


(A.  4a) 
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and 


U 


d* 


(A.  4b) 


(A.  5a) 


These  terms  may  be  substituted,  as  appropriate,  in  (A.l)  to  obtain: 
where  the  continuity  equation  is  given  by, 


(A.  5b) 


(A.  6) 


It  should  be  noted  that  (A.  6)  differs  slightly  from  the  usual  flux  form 
in  which  the  map  factor  appears  in  both  flux  quantities.    Thus,  the 
initial  divergence  will  not  vanish  identically  (after  use  of  the  balance 
equation) , 
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APPENDIX  B 
Energy  Computations 

In  pressure  coordinates  the  total  potential  energy  (internal  plus 
geopotential)  may  be  written   (for  a  smooth  earth)  as, 

A  o 

where  "dA"  represents  an  area  increment.     Now,  if  one  eliminates 
from  the  hydrostatic  equations  in  pressure  and  zeta  coordinates, 
respectively,  one  derives  the  relationship, 

Si    =     -    -£.<??    =     -2   JC  (B.2) 

r  s<r  S 

where  "s"  is  the  vertical  mapping  factor,  and  1T  is  the  lower  boundary 
pressure  (sea  level).    Substituting  (B.2)  in  the  first  expression,  the 
total  potential  energy  becomes, 

r  _   Sif  (   irr   ^j  dA  (b.3) 

p      3   J  )     s 

A  o 

In  a  similar  manner,  the  appropriate  form  for  kinetic  energy  becomes, 

Eve   t%  tf|(uS^  drdA  (B"4) 

Ao 
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The  corresponding  finite -difference  forms  for  (B  ,  3)  and  (B.4)  were 
given  as , 

and 


(B.6) 


where  "N"  is  the  number  of  mesh  areas  and  "L"  is  the  number  of 
layers.    In  addition,   "D"  is  the  grid  mesh  length,  and  other  symbols 
were  as  initially  specified. 
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APPENDIX    C 
Divergence  Equation  Forms  in  the  Sigma  System 

One  can  take  the  divergence  of  the  vector  equation  of  frictionless , 
adiabatic  motion  on  a  sigma  surface  and  show  that  the  resulting 
divergence  equation  assumes  the  form, 

4  +  §l  -*£  -  sa^o  +  vo-.^  +  v$  -u*£ 

d-fc  f  dc  $x 


**$  *  ^-wn^H™^ 


=  0(C.l) 


if  the  "sphericity"  terms  are  omitted  from  the  equation  of  motion. 

It  has  been  shown  (see  Haltiner  (3),  for  example)  that  the  terms  of 
significance  to  less -generalized  vorticity  models  specialize  to 


v^=  -*-iv§-^v+[ 


(C.2) 


usually  referred  to  as  the  linear  balance  equation,   under  the  assumption 
that  the  non-divergent  component  of  the  wind  is  large  compared  to  the 
divergent  component. 

In  sigma  coordinates,  however,  terms  arising  from  the  additional 
pressure  force  term  in  the  momentum  equations  are  rather  large,  and 
must  be  included  in  the  balance  equation  in  sigma  coordinates.    The 
resulting  form  is  given  as 


"3x J  +  I 
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=  o        (°-3) 


If  one  applies  the  (lk*  Vx      )  operator  to  the  same  equation  of 
motion,  the  vertical  component  of  the  relative  vorticity  on  a  sigma 
surface  is  found  to  be 


^i^+Knv;'"i^)+ff), 


T 


(C.4) 


which  preserves  the  mathematical  form  of  the  linear  balance  equation. 
It  should  be  noted  that  (C .  3)  takes  on  a  more  complicated  form  when 
the  term 

"T —-  (u,v)  1? 

is  added  to  the  momentum  equations  . 


*    as  appropriate 


59 


THESIS  DISTRIBUTION  LIST 

Copies 

1.  LT  Philip  G.  Kesel,  USN  5 
Fleet  Numerical  Weather  Central 

NPGS,   Monterey,  California  93940 

2.  Professor  G.J.  Haltiner  3 
Department  of  Meteorology 

NPGS,   Monterey,  California  93940 

3.  Professor  R.T.  Williams  3 
Department  of  Meteorology 

NPGS,  Monterey,  California  93940 

4.  Library,  Naval  Postgraduate  School  2 
Monterey,  California  93940 

5.  Defense  Documentation  Center  20 
Cameron  Station 

Alexandria,  Virginia  22314 

"'    6.         Naval  Weather  Service  Command  1 

Washington  Navy  Yard 
Wasington,  D.C.   20390 

^7.  Officer  in  Charge,  NWRF  1 

Naval  Air  Station,  Bldg.  R-48 
Norfolk,  Virginia  23511 

8.  Commanding  Officer  2 
FNWC,  NPGS 

Monterey,  California  93940 

9.  AFCRL  -  Research  Library  1 
L.G.   Hanscom  Field  (Attn:  Nancy  Davis,  Stop  29) 
Bedford,   Massachusetts  01730 

10.  Dr.  G.B.  Tucker  1 
Commonwealth  Bureau  of  Meteorology 

P.O.   Box  1289  K 
Melbourne,  Vic,  Australia 

11.  Department  of  Meteorology  (Code  51)  3 
Naval  Postgraduate  School 

Monterey,  California  93940 


60 


12.  Department  of  Oceanography  (Code  58) 
Naval  Postgraduate  School 
Monterey,  California  93940 

13.  American    Meteorological  Society 
4  5  Beacon  Street 

Boston,   Massachusetts  02128 

14.  Commander,  Air  Weather  Service  (MAC) 
U.S.  Air  Force 

Scott  AFB,  Illinois  62226 

15.  Atmospheric  Sciences  Library 
Environmental  Science  Service  Administration 
Silver  Spring,  Maryland  20910 

16.  Oceanographer  of  the  Navy 
The  Madison  Building 

732  N.  Washington  Street 
Alexandria,  Virginia  22314 

'    17.        Naval  Oceanographic  Office 
Attn:     Library 
Washington,  D.C.   20390 

18.  Naval  Oceanographic  Data  Center 
Washington,  D.C.   20390 

19.  Director,  Maury  Center  for  Ocean  Sciences 
Naval  Research  Laboratory 
Washington,  D.C.   20390 

20.  Dr.  A.  Arakawa 

UCLA  Department  of  Meteorology 
Los  Angeles,  California  90024 

21.  Professor  N.A.   Phillips 
MIT  54-1422 

Cambridge,  Massachusetts  02139 

22.  Dr.  F.H.   Bushby 
Meteorology  Office 
Bracknell,  England 

23.  Dr.  F.G.  Shuman 
Director,  NMC  ,  ESSA 
Suitland,  Maryland 


61 


UNCLASSIFIED 


Security  Classification 


DOCUMENT  CONTROL  DATA  •  R&D 

(Security  claaaltlcatlon  of  till*,   body  of  mbtumct  and  Indaxfng  annotation  muat  bo  anitrad  whan  tfia  oratatt  report  la  claaaillad) 


1      ORIGIN  A  TIN  C  ACTIVITY  (Corporata  author) 

Naval  Postgraduate  School 
Monterey,  California  93940 


2a.    KI^OKT    IICUNITY     C    L  ASS  I  F  I  C  A  T  I  O  r 


Unclassified 


2  6    anou* 


3     REPORT  TITLE 

EXPERIMENTS  WITH  ATMOSPHERIC  PRIMITIVE -EQUATION  MODELS 


4     DESCRIPTIVE  NOTES  (Typa  of  report  and  fnelu 

Master's  Thesis,  1968 


S     AUTHORfS;  (Laat  nam.     Om  mni,    initial) 

KESEL,   Philip  G.,  LT,   USN 


•    REPO  RT  DATE 

December  1968 


7«      TOTAL  NO.   OF    >*tll 


60 


7*.    NO.   OF  RIPS 


20 


8a      CONTRACT    OR    ORANT    NO. 


6.    PROJ1CT   NO. 


•  a.    ORIOINATOR't   RtRORT   NUMSIW'Jj 


tb.   QTHKR  RIRORT    no(S)   (A  ny  othat  numbars  fl«l  may  ba  aaalgnad 
mia  taport) 


10    A  V  A  I  LABILITY/LIMITATION  NOTICES 

Distribution  of  this  document  is  unlimited. 


11-  SUPPLEMENTARY  NOTES 


12    SPONSORING  MILITARY  ACTIVITY 

Naval  Postgraduate  School 
Monterey,  California  93940 


13     ABSTRACT 

Two  atmospheric  prediction  models,  based  upon  the  meteorological 
primitive  equations,  were  programmed  and  tested  on  the  CDC  6500  computer. 
A  frictionless ,  barotropic  model  was  integrated,  using  ten-minute  time  steps, 
at  500  MBS  for  periods  up  to  four  days.    Then,  an  n-level,  baroclinic  model 
for  an  inviscid,  adiabatic  atmosphere  was  integrated  for  brief  test  periods  at 
five  uniformly-spaced  zeta  (sigma)  surfaces  over  a  smooth  earth.    Both  models 
were  initialized  with  actual  data  fields  provided  by  FNWC  Monterey. 

Acceptability  was  based  on  measurements  of  energy  and  vorticity 
parameters,  as  well  as  on  qualitative  assessments  of  output  fields.    The 
barotropic  500  MB  height  forecasts  were  found  to  be  comparable  to  the  FNWC 
(vorticity  model)  barotropic  forecasts.    Mean-square-vorticity  and  kinetic 
energy  were  suitably  conserved  for  integrations  up  to  four  days.    An  energy 
accumulation  in  the  smaller  range  of  scale  was  increasingly  evident  beyond 
two  days  into  the  forecast  period.    No  smoothing  was  performed  to  control 
this  accumulation.    A  limited  number  of  test  integrations  was  made  with 
the  five -level  model.    Computational  instabilities  were  observed  for  forecasts 
beyond  twelve  hours.    Interim  results  are  presented. 


DD 


FORM 

1    JAN    64 


1473 


UNCLASSIFIED 


63 


Security  Classification 


UNCLASSIFIED 


Security  Classification 


KEY     WORDS 


PRIMITIVE-EQUATION  MODELS 
BAROTROPIC  P.E.    MODELS 
BAROCLINIC  P.E.   MODELS 
INTEGRATION  OF  ATMOSPHERIC  PRIMITIVE- 
EQUATIONS 


LINK     B 


DD,?0R:.,1473  (back) 


UNCLASSIFIED 


/•.        101  •■07.1 


64 


Security  Classification 


f/otii  — — 


Ljaulorc 

_JSHEIF   BINDER 

^ZZT      Syracuse,  N.  Y. 
|      Slocklon,  Calif. 


thesK393 


3  2768  00416533  2 

^^HJDLEY  KNOX  LIBRARY 


